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Abstract Financial derivatives pricing aims to find the fair value of a financial 
contract on an underlying asset. Here we consider option pricing in the partial dif¬ 
ferential equations framework. The contemporary models lead to one-dimensional 
or multidimensional parabolic problems of the convection-diffusion type and gen¬ 
eralizations thereof. An overview of various operator splitting methods is presented 
for the efficient numerical solution of these problems. 

Splitting schemes of the Alternating Direction Implicit (ADI) type are discussed 
for multidimensional problems, e.g. given by stochastic volatility (SV) models. For 
jump models Implicit-Explicit (IMEX) methods are considered which efficiently 
treat the nonlocal jump operator. Eor American options an easy-to-implement oper¬ 
ator splitting method is described for the resulting linear complementarity problems. 

Numerical experiments are presented to illustrate the actual stability and conver¬ 
gence of the splitting schemes. Here European and American put options are consid¬ 
ered under four asset price models: the classical Black-Scholes model, the Merton 
jump-diffusion model, the Heston SV model, and the Bates SV model with jumps. 
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1 Introduction 

In the contemporary international financial markets option products are widely 
traded. The average daily turnover in the global over-the-counter derivatives mar¬ 
kets is huge. For example, in the foreign exchange market this was approximately 
equal to 337 billion US dollars in April 2013 [5]. In addition to standard call and 
put options, the so-called vanilla options, a broad range of exotic derivatives ex¬ 
ists. One of the primary goals of financial mathematics is to determine the fair val¬ 
ues of these derivatives as well as their sensitivities to underlying variables and 
parameters, which are crucial for hedging. To this purpose, advanced mathemati¬ 
cal models are employed nowadays, yielding initial-boundary value problems for 
time-dependent partial differential equations (PDFs) and generalizations thereof, 
see e.g. [4, 14, 59, 75, 77, 85]. These problems are in general multidimensional and 
of the convection-diffusion kind. In some cases analytical formulas in semi-closed 
form for the exact solutions have been obtained in the literature. For the majority 
of option valuation problems, however, such formulas are not available. In view of 
this, one resorts to numerical methods for their approximate solution. To banks and 
other financial institutions, the efficient, stable, and robust numerical approximation 
of option values and their sensitivities is of paramount importance. 

A well-known and versatile approach to the numerical solution of time-dependent 
convection-diffusion equations is given by the method of lines. It consists of two 
general, consecutive steps. In the first step the PDE is discretized in the spatial vari¬ 
ables, e.g. by finite difference, finite volume, or finite element methods. This leads 
to a so-called semidiscrete system of ordinary differential equations. In the second 
step the obtained semidiscrete system is numerically solved by applying a suitable, 
implicit time-discretization method. If the PDE is multidimensional, then the latter 
task can be computationally very intensive when standard application of classical 
implicit methods, such as the Crank-Nicolson scheme, is used. In the recent years, 
a variety of operator splitting methods have been developed that enable a highly ef¬ 
ficient and stable numerical solution of semidiscretized multidimensional PDEs and 
generalizations thereof that arise in financial mathematics. 

The aim of this chapter to give an overview of main classes of operator splitting 
methods with applications in finance. Here we have chosen to consider a variety 
of, increasingly sophisticated, models that are well-known in the financial option 
valuation literature. 

We deal in the following with two basic types of options, involving a given so- 
called strike price K > 0 and a given maturity time T > 0, where today is always 
denoted by time 0. A European call (put) option is a contract between two parties, 
the holder and the writer, which gives the holder the right to buy from (sell to) the 
writer a prescribed asset for the price K at the future date T. An American call (put) 
option is the same, except that the holder can exercise at any time between today and 
the maturity date. An option is a right and not an obligation. The underlying asset 
can be a stock, a foreign currency, a commodity, etc. For a detailed introduction to 
financial options we refer to [45]. Clearly, an option has value and a central question 
in financial mathematics is what its fair value is. 
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2 Models for Underlying Assets 
2.1 Geometric Brownian Motion 

The seminal papers by Black & Scholes [7] and Merton [63] present a key equation 
for the fair values of European call and put options. In these papers the dynamics of 
the underlying asset price is modeled by the stochastic differential equation (SDE) 


dS(t) = ^S{t)dt + <jS{t)dW{t) (f>0). 


( 1 ) 


Here W{t) denotes the Wiener process or standard Brownian motion, and /r, cj are 
given real parameters that are called the drift and the volatility, respectively. The 
volatility is a degree for the uncertainty of the return realized on the asset. 

The SDE (1) describes a so-called geometric Brownian motion, which satisfies 
5(f) > 0 whenever 5(0) > 0. Under this asset price model and several additional as¬ 
sumptions, Black, Scholes, and Merton derived the famous partial differential equa¬ 
tion (PDE) 



( 2 ) 


Here u{s,t) represents the fair value at time T — f of a European vanilla option if 
S{T — t) = s. The quantity r in (2) is the risk-free interest rate and is given. A main 
consequence of the Black, Scholes, and Merton analysis is that the drift ji actually 
does not appear in the option pricing PDE. This observation has led to the important 
risk-neutral valuation theory. It is beyond the scope of the present chapter to discuss 
this theory in more detail, but see e.g. [45, 75]. 

In formulating (2) we have chosen f as the time till maturity. Thus the time runs 
in the opposite direction compared to (1). Accordingly, the payoff function 0, which 
defines the value of the option contract at maturity time T, leads to an initial condi¬ 
tion 


m(s,O) = 0(s) (s>0). 

Eor a European vanilla option with given strike price K there holds 


(3) 



max(s —Ai,0) for s>0 (call), 
max(A' —s,0) for s>0 (put). 


(4) 


and at i = 0 one has the Dirichlet boundary condition 



for 0<t<T (call), 
for 0 < f < T (put). 


(5) 


Equation (2) is called the Black-Scholes PDE or Black-Scholes-Merton PDE. It 
is fully deterministic and it can be viewed as a time-dependent convection-diffusion- 
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reaction equation. For European vanilla options, an analytical solution u in semi- 
closed form was derived in [7], constituting the well-known Black-Scholes formula. 

The Black-Scholes PDE is generic in the sense that it is valid for a wide range of 
European-style options. The initial and boundary conditions are determined by the 
specific option. As an example, for a European up-and-out call option with given 
barrier B > K, the PDE (2) holds whenever 0 < s < B, 0 < f < T. In this case, the 
initial condition is 


u{s,0) =max{s — K,0) for 0 < i < B 
and one has the Dirichlet boundary conditions 

u{Q,t) = u{B,t) = 0 for 0 <t <T. 

The homogeneous condition at i = B corresponds to the fact that, by construction, 
an up-and-out call option becomes worthless whenever the underlying asset price 
moves above the barrier. 

Eor many types of options, including (continuous) barrier options, semi-analytical 
pricing formulas have been obtained in the literature in the Black-Scholes frame¬ 
work, see e.g. [45]. At present it is well-known, however, that each of the assump¬ 
tions underlying this framework are violated to a smaller or larger extent in practice. 
In particular, the interest rate r and the volatility <7 are not constant, but vary in time. 
In view of this, more advanced asset pricing models have been developed and, as 
a consequence, more advanced option valuation PDEs are obtained. In this chap¬ 
ter we do not enter into the details of the mathematical connection between asset 
price SDEs and option valuation PDEs, but mention that a main tool is the cele¬ 
brated Eeynman-Kac theorem, see e.g. [75]. In the following we discuss typical, 
contemporary instances of more advanced option valuation PDEs. 


2.2 Stochastic Volatility and Stochastic Interest Rate Models 


Heston [38] modeled the volatility itself by a SDE. The Heston stochastic volatility 
model is popular especially in the foreign exchange markets. The corresponding 
option valuation PDE is 


du I 9 d^u d^u I 9 d^u du , 

— = + pasv^:-^ + ka + rs— y k{7] 

at ds^ dsdv dv ds 



— rii 


( 6 ) 


for i > 0, V' > 0, and 0 <t <T. Here u{s, v', t) represents the fair value of a European- 
style option if at t time units before maturity the asset price equals s and the vari¬ 
ance equals v. We note that by definition the variance is the square of the volatility. 
The positive parameters K and rj are the mean-reversion rate and long-term mean, 
respectively, of the variance. O’ > 0 is the volatility-of-variance, and p € [—1,1] 
denotes the correlation between the two underlying Brownian motions. Equation 
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(6) is called the Heston PDE. It can be viewed as a time-dependent convection- 
diffusion-reaction equation on an unbounded, two-dimensional spatial domain. If 
the correlation p is nonzero, which almost always holds in practice, then the Heston 
PDE contains a mixed spatial derivative term. 

For a European vanilla option under the Heston model, one has an initial con¬ 
dition as well as a boundary condition at s = 0 that are the same as in the Black- 
Scholes case discussed above. In the Heston case there is also a boundary v = 0. 
Observe that as v | 0, then all second-order derivative terms vanish in (6). It has 
been proved in [25] that for the fair option value function u the Heston PDE is 
fulfilled if V' = 0, which constitutes the (nonstandard) boundary condition at v' = 0. 

For the Heston asset pricing model (which we did not explicitly formulate) the 
so-called Feller condition 2Kr\>a^ is often considered in the literature. This condi¬ 
tion determines whether or not the variance process can attain the value zero (given 
a strictly positive initial variance): it cannot attain zero if and only if Feller holds. 
The situation where the Feller condition is violated is well-known to be challenging 
when numerically solving the Heston asset pricing model. For the Heston option 
valuation PDE (6), on the other hand, it turns out that this issue is not critical in the 
numerical solution. 

A refinement of the Heston model is obtained by considering also a stochastic 
interest rate, see e.g. [32, 33, 35, 36]. As an illustration we consider the case where 
the interest rate is described by the well-known Hull-White model [45, 46]. This 
leads to the following so-called Heston-Hull-White (HHW) PDE for the option 
value function u = u{s,v,r,t): 


dll 


2" + 2^1 + 2‘^2 ^ +P12(T1«'^ + Px^a2S^V — 

_ d^ii 

-FP23CTi(J2VV' 


dvdr 


dll , ^ du , , , , ( 9 m 

rs— 4- K(ri -v)—+a(b(T-t)-r)- - ru 

ds dv dr 


(7) 


for i > 0, V' > 0, —oo < r < °o, and 0 < t <T. Here K, rj, <Ji, a, and 02 are given 
positive real constants and b denotes a given deterministic, positive function of 
time. Further, there are given correlations pi 2 , P 13 , P 23 G [—1, !]■ Clearly, the HHW 
PDE is a time-dependent convection-diffusion-reaction equation on an unbounded, 
three-dimensional spatial domain with three mixed derivative terms. For a Euro¬ 
pean vanilla option, initial and boundary conditions are the same as in the Heston 
case above. Note that if v 1 0, then all second-order derivative terms, apart from the 
d^u/dr^ term, vanish in (7). 

The Heston and HHW models are two of many instances of asset pricing models 
that lead to multidimensional option valuation PDEs. Multidimensional PDFs are 
also obtained when considering other types of options, e.g. options on a basket of 
assets. Then, in the Black-Scholes framework, the dimension of the PDE is equal to 
the number of assets. In general, analytical solutions in (semi-)closed form to these 
PDEs are not available. 
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2.3 Jump Models 


Sometimes the value of the underlying asset changes so rapidly that this would have 
very tiny probability under the above Brownian motion based models. For example, 
the stock price during a market crash or after a major news event can move very fast. 
Already in 1976, Merton proposed in [64] to add a jump component in the model 
of the underlying asset price. In his model, the jumps are log-normally distributed 
and their arrival times follow a Poisson process. After a jump the value of the asset 
is obtained by multiplying the value before the jump by a random variable with the 
probability density function (PDF) 


m 


1 


ydV^ 


exp 


(logy- 7 ) 2 \ 

252 J 


( 8 ) 


for y > 0 , where 7 is the mean of the normal distribution and 5 is its standard de¬ 
viation. Kou proposed in [56] a log-double-exponential distribution defined by the 


PDF 


fiy) 


^a27“^ 0 <y<l, 

paiy-“'-', y>l, 


(9) 


where p,q,a\ > 1, and ai are positive constants such that p + q = 1. These models 
have finite jump activity which is denoted by X here. There are also many popular 
infinite jump activity models like the CGMY model [11]. In the following we shall 
consider only finite activity models. 

The value u{s,t) of a European option satisfies the partial integro-differential 
equation (PIDE) 


— = \a^s^-^X{r-XQs—-{r + X)u + Xj^ u{sy,t)f{y)dy (10) 
for s > 0 and 0 < t < T, where ^ is the mean jump size given by 

C= r{y-my)dy. (H) 

Jo 

Eor the Merton and Kou models the mean jumps are = e1'+^^/2 — 1 and + 

— 1 , respectively. 

Bates proposed to combine the Heston stochastic volatility model and the Merton 
jump model in [ 6 ]. Under this model the value u{s,v,t) of a European option satisfies 
the PIDE 


du 


= ii2y 


d^u 


pasv 


d^u 

dsdv 


1 n d U , . y, y u It . 


du 

'Ts 


-{r + X)u + X u{syyV,t)f{y)dy 

Jo 



(12) 
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for s > 0, V > 0, and 0 < t < T, where the PDF / is given by (8). For an extensive 
discussion on jump models in finance see e.g. [16]. 


3 Linear Complementarity Problem for American Options 

Unlike European-style options, American-style options can be exercised at any time 
up to the maturity date. Hence, the fair value of an American option is always greater 
than or equal to the instantaneous payoff. 


u> (j). 


(13) 


Due to this early exercise constraint, the P(I)DE does not hold everywhere anymore. 
Instead, a linear complementarity problem (LCP) or partial (integro-)differential 
complementarity problem is obtained in general for the fair value of an American 
option: 



(14) 


where £/ stands for the pertinent spatial differential operator. Eor example, for the 
Black-Scholes model. 


1 2 2 

MM — —<y s + rs^ - ru. 

2 os^ os 


The above inequalities and equation hold pointwise. The equation in (14) is the 
complementarity condition. It states that at each point one of the two inequalities has 
to be an equality. The paper [44] discusses the LCP formulation for American-style 
options under various asset price models and studies the structure and properties of 
the obtained fully discrete LCPs. 

We note that the penalty approach is a popular alternative for LCPs. Here a 
penalty term is added to the P(I)DE for a European option with the aim to enforce 
the early exercise constraint (13). The resulting problems are nonlinear and their 
efficient numerical solution is considered in [27], for example. Eor several other 
alternative formulations and approximations for LCPs, we refer to [80]. 


4 Spatial Discretization 


In this chapter we employ finite difference (ED) discretizations for the spatial deriva¬ 
tives. An alternative approach would be to use finite element discretizations; see e.g. 
[1, 74]. It is common practice to first truncate the infinite s-domain [0,°°) to [Oj^max] 
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with a sufficiently large, real 5max- Typically one wishes iSmax to be such that the er¬ 
ror caused by this truncation is a small fraction of the error due to the discretization 
of the differential (and integral) operators. Similarly, with multidimensional models 
including the variance v or the interest rate r, their corresponding infinite domains 
are truncated to sufficiently large bounded domains. The truncation requires addi¬ 
tional boundary conditions to be specified. For an actual choice of these conditions 
for the models considered in Sections 2, 3 we refer to Section 7. 

Let the grid in the s-direction be defined by the mi -f 1 grid points 0 = so < 
Si < • • • < Sm, = 5max- The Corresponding grid sizes are denoted by Asi = Si — s,_i, 
i= 1,2,..., mi. For multidimensional models, we use tensor product grids. For ex¬ 
ample, in the case of a stochastic volatility model, if a grid for the variance v is given 
by 0 = I'o < I'l < • • • < v,n 2 = Fmax, then (mi 4- 1) x (m 2 + 1) spatial grid points are 
defined by (s,, v/) with / = 0,1,... ,mi and y = 0,1,... ,m 2 . In financial applica¬ 
tions nonuniform grids are often preferable over uniform grids. The use of suitable 
nonuniform grids will be illustrated in Section 7. 

For discretizing the first derivative ^ and the second derivative at s = s,, 
we employ in this chapter the well-known central FD schemes 


dui 


-4i,+i 


4s;(4si -|-4s,+ l) 


Asi+i-Asi 

-1 H-W i 

AsiAsi+i 


Asi 


(4s;-|-4si+i)4si+i 




(15) 


and 


d^Ui 2 

ds^ Asi{Asi + Asi+i) 


Ui-l - 


ZiiiZii,'+l (zii,'-|-ziii+i)zii/+i 


Ui+i- (16) 


With multidimensional models the analogous schemes are used for the other spatial 


dll j 


dhii 


d~U' 

at V = Vj. For the mixed derivative 


at 


directions, thus e.g. for and 
(i,v) = {si,Vj) we consider the 9-point stencil obtained by successively applying 
the central FD schemes for the first derivative in the s- and v-directions. With suffi¬ 
ciently smooth varying grid sizes, the above central FDs give second-order accurate 
approximations for the derivatives. 

We mention that in financial applications other FD schemes are employed as 
well, such as upwind discretization for first derivative terms or alternative discretiza¬ 
tions for mixed derivative terms. 

With the jump models the integral term needs to be discretized at grid points s, . 
First the integral is divided into two parts 

1'°° rSmax/^i r°° 

/ u{siy,t)f{y)dy= i u{siy,t)f{y)dy+ u{siy,t)f{y)dy, 

JQ Jo JSmax/^i 


which correspond to the values of u in the computational domain [Oj^max] and out¬ 
side of it, respectively. The second part can be estimated using knowledge about u 
in the far field [5'max,°°)- For example, for put options u is usually assumed to be 
close to zero for s > 5max and, thus, the second integral is approximated by zero in 
this case. The PDFs / are smooth functions apart from the potential jump at y = 1 
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in the Kou model. Due to the smoothness of the integrand the trapezoidal rule leads 
to second-order accuracy with respect to the grid size. This gives the approximation 


f SmSix/^i 

/ u{siy,t)f{y)dy', 

JQ 


"'1 A ?- 
;=i 




For example, the papers [71] and [78] describe more accurate quadrature rules for 
the Merton and Kou jumps models, respectively. The discretization of the integral 
term leads to a dense matrix. The integral can be transformed into a convolution 
integral and due to this FFT can be used to compute it more efficiently; see [2, 3, 
22, 77], for example. In the case of the Kou model, efficient recursion formulas can 
be used [12, 78]. 


5 Time Discretization 
5.1 The 9-method 

For any P(I)DE from Section 2, the spatial discretization outlined in Section 4 leads 
to an initial value problem for a system of ordinary differential equations, 

t/(f)=A(f)[/(f) + G(f) {Q<t<T), U{Q) = Uo. (17) 

Here A(f) for 0 < f < T is a given square real matrix and G{t) is a given real vector 
that depends on the boundary conditions. The entries of the solution vector U{t) 
represent approximations to the exact solution of the option valuation P(I)DE at the 
spatial grid points, ordered in a convenient way. The vector Uq is given by direct 
evaluation of the option’s payoff function at these grid points. 

The semidiscrete system (17) is stiff in general and, hence, implicit time dis¬ 
cretization methods are natural candidates for its numerical solution. Let parameter 

6 € (0,1] be given. Let time step At = T/N with integer A > 1 and temporal grid 
points t„ = nAt for integers Q<n<N. The 6-method forms a well-known implicit 
time discretization method. It generates approximations U„ to U (f„) successively for 
n= 1 , 2 ,... ,A by 


Un = Un-x + (1 - 0)Af A(f„_i)t /„_1 + 0Af A(f„)t/„ + Af G„_i+e, (18) 

where G„_i +0 denotes an approximation to G(f) at f = (n — 1 -F 6)At. This can also 
be written as 


(I - 0AfA(G))G„ = (!+(!- 0)Af A(f„_i))G„_i +Af G„_i+ 0 , 

with I the identity matrix of the same size as A(f). Eor 0 = 1 one obtains the first- 
order backward Euler method and for 0 = 5 the second-order Crank-Nicolson 
method or trapezoidal rule. Eor simplicity we consider in this chapter only con- 
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Slant time steps, but most of the presented time discretization methods can directly 
be extended to variable time steps. 

When applying the Crank-Nicolson method, it is common practice in finance to 
first perform a few backward Euler steps to start the time stepping. This is often 
called Rannacher smoothing [67]. It helps to damp high-frequency components in 
the numerical solution, due to the nonsmooth initial (payoff) function, which are 
usually not sufficiently damped by the Crank-Nicolson method itself. 

Clearly, in order to compute the vector t/„ defined by (18), one has to solve a 
linear system of equations with the matrix I— 9AtA{t„). When the option valuation 
PDE is multidimensional, the size of this matrix is usually very large and it possesses 
a large bandwidth. Eor a PIDE, this matrix is dense. In these situations, the solution 
of the linear system can be computationally demanding when standard methods, 
like LU decomposition, are applied. Time discretization methods based on operator 
splitting can then form an attractive alternative. The key idea is to split the matrix 
A(f) into several parts, each of which is numerically handled more easily than the 
complete matrix itself. 


5.2 Operator Splitting Methods Based on Direction 

Eor multidimensional PDEs, splitting schemes of the Alternating Direction Implicit 
(ADI) type are often applied in financial practice. To illustrate the idea, the two- 
dimensional Heston PDE and three-dimensional HHW PDE, given in Section 2.2, 
are considered. Eor the Heston PDE the semidiscrete system (17) is autonomous; 
we split 

A = Aq -t- Ai -f A 2 . 

Next, for the HHW PDE, 


A(f) = Aq + Ai -f A2-I-A3(r). 

Here Aq is chosen as the part that represents all mixed derivative terms. It is nonzero 
whenever (one of) the correlation factor(s) is nonzero. The parts Ai, A 2 , and A 3 (f) 
represent all spatial derivatives in the s-, v-, and r-directions, respectively. The latter 
three matrices have, possibly up to permutation, all a fixed small bandwidth. The 
vector G(f) in the semidiscrete system is splitted in a similar way. Eor notational 
convenience, define functions F, by 

F,(f,y)=A,y + G,- (; = 0,1,2) and F 3 (f,y) = A 3 (f)y+ G 3 (f) 

for 0 < f < r, y e R"”. Set F = 1^=0 Fy with k = 2 for Heston and k = 3 for HHW. 
We discuss in this section four contemporary ADI-type splitting schemes; 
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Douglas (Do) scheme 

io = Urt-l 

< Yj = Yj^i + eAtiFj{tn,Yj)-Fj{t„^uUn-i)) (19) 

^ Un = Yk. 

Craig-Sneyd (CS) scheme 

Yj = Yj^i + 9At (Fj{tn,Yj)- Fjitn^ i, t/„-1)) (; = 1,2,...,, 

< Fq = lo+5^1 (Fo(f„,Fi:) — Fo(f„_i,t/„_i)), (20) 

Yj = Yj^i + eAt (Fj (tn ,Yj)- Fjitn^ 1, t/„-1)) (j = 1,2,..., k), 

. Un = %. 

Modified Craig-Sneyd (MCS) scheme 
' Yq = t/„_i 

Yj = Yj^i + eAtiFjitn,Yj)-Fjit„^i,Un-i)) (j =1,2,.. .,k), 

Yo = Yo+eAtiFoitn,Yk)-Fo{tn-l,Un-l)), 

< ^ (21) 
Fo = Fo + (i - 0)4f (F (tn , Fi) - F(f„_ i, t/„-1)), 

Yj = Yj^i + 9At (Fj(tn,Yj) — Fy(f„_i, t/„_i)) (j = i,2,... ,k), 

. U„ = %. 

Hundsdorfer-Verwer (HV) scheme 

' Fq = t/„_l + AtF{tn-l,Un-l), 

Yj = Yj^i + 9At (Fj(tn,Yj)- Fj(tn- 1, t/„-1)) (; = 1,2,...,, 

< ?o = Fo + i4f(F(f„,y^)-F(f„_i,t/„_i)), (22) 

Yj = Yj^, + 9At(Fj(tn,Yj)-Fj(tn,Yk)) (j = \,2,... ,k), 

. Un = %. 
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In the Do scheme (19), a forward Euler predictor step is followed by k implicit 
but unidirectional corrector steps that serve to stabilize the predictor step. The CS 
scheme (20), the MCS scheme (21), and the HV scheme (22) can be viewed as 
different extensions to the Do scheme. Indeed, their first two lines are identical to 
those of the Do scheme. They next all perform a second predictor step, followed by k 
unidirectional corrector steps. Observe that the CS and MCS schemes are equivalent 
if (and only if) 0 = 5 . 

Clearly, in all four ADI schemes the Ao part, representing all mixed derivatives, 
is always treated in an explicit fashion. In the original formulation of ADI schemes 
mixed derivative terms were not considered. It is a common and natural use in the 
literature to refer to the above, extended schemes also as ADI schemes. In the special 
case where Fq = 0, the CS scheme reduces to the Do scheme, but the MCS scheme 
(with 9 ^ j) and the HV scheme do not. Following the original ADI approach, the 
Ai, A 2 , A 3 (f) parts are treated in an implicit fashion. In every step of each scheme, 
systems of linear equations need to be solved involving the matrices (I — 9 At A j) 
for j = 1,2 as well as (I — 9At A^{t„)) if A: = 3. Since all these matrices have a 
fixed, small bandwidth, this can be done very efficiently by means of LU decom¬ 
position, cf. also Section 6.1. Because for j = 1,2 the pertinent matrices are further 
independent of the step index 11 , their LU decompositions can be computed once, 
beforehand, and then used in all time steps. Accordingly, for each ADI scheme, 
the number of floating point operations per time step is directly proportional to the 
number of spatial grid points, which is a highly favorable property. 

By Taylor expansion one obtains (after some elaborate calculations) the classical 
order of consistency* of each ADI scheme. For any given 0, the order of the Do 
scheme is just one whenever Aq is nonzero. This low order is due to the fact that 
the Ao part is treated in a simple, forward Euler fashion. The CS scheme has order 
two provided 0 = 5 . The MCS and HV schemes are of order two for any given 0. 
A virtue of ADI schemes, compared to other operator splitting schemes based on 
direction, is that the internal vectors Yj, Yj form consistent approximations to U{t„). 

The Do scheme can be regarded as a generalization of the original ADI schemes 
for two-dimensional diffusion equations by Douglas & Rachford [23] and Peaceman 
& Rachford [ 66 ] to the situation where mixed derivative terms are present. This gen¬ 
eralization was first considered by McKee & Mitchell [61 ] for diffusion equations 
and subsequently in [62] for convection-diffusion equations. 

The CS scheme was developed by Craig & Sneyd [18] with the aim to obtain 
a stable second-order ADI scheme for diffusion equations with mixed derivative 
terms. 

The MCS scheme was constructed by In’t Hout & Welfert [43] so as to arrive at 
more freedom in the choice of 0 as compared to the second-order CS scheme. 

The HV scheme was designed by Hundsdorfer [47] and Verwer et. al. [83] for 
the numerical solution of convection-diffusion-reaction equations arising in atmo¬ 
spheric chemistry, cf. also [48]. The application of the HV scheme to equations 
containing mixed derivative terms was first studied in [42, 43]. 


* That is, the order for fixed nonstiff ODE systems. 
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The Do and CS schemes are well-known for PDEs in finance, see e.g. [4, 59]. 
More recently, the MCS and HV schemes have gained interest, see e.g. [14, 20, 24, 
35, 36, 39, 54]. 

The formulation of the ADI schemes (19)-(22) is analogous to the type of for¬ 
mulation used in [47]. In the literature, ADI schemes are also sometimes referred to 
as Stabilizing Correction schemes, and are further closely related to Approximate 
Matrix Factorization methods and Implicit-Explicit (IMEX) Runge-Kutta methods, 
cf. e.g. [48]. 

In [40,41,42,43] comprehensive stability results in the von Neumann sense have 
been derived for the four schemes (I9)-(22) in the application to multidimensional 
convection-diffusion equations with mixed derivative terms. These results concern 
unconditional stability, that is, without any restriction on the time step At. For each 
ADI scheme, lower bounds on 9 guaranteeing unconditional stability have been ob¬ 
tained, depending in particular on the spatial dimension. Based on these theoretical 
stability results and the numerical experience in [35, 36, 39] the following values 
are found to be useful for k = 2,3: 

• Do scheme with 6 = ^ (if k = 2) and 0 = | (if k = 3) 

• CS scheme with 0 = j 

• MCS scheme with 0 = ^ (if k = 2) and 0 = maxl^, ■^{2y+ 1)} (if k = 3) 

• HV scheme with 0 = i + ^V3. 

Here 7 = max{|pi 2 |, IPbI, IP23I} C [0,1], which is a measure for the relative size of 
the mixed derivative coefficients. 

In addition to ADI schemes, there exists a variety of well-known alternative 
operator splitting schemes based on direction, called Locally One-Dimensional 
(LOD) methods, fractional step methods, or componentwise splitting schemes. 
These schemes originate in the 1960s in the work by Dyakonov, Marchuk, Samarskii, 
Yanenko, and others. Some of them are related to Strang splitting schemes, devel¬ 
oped at the same time. For a general overview and analysis of such methods we refer 
to [48, 60]. Applications in financial mathematics of these schemes are considered 
in, for example, [50, 79]. 


5.3 Operator Splitting Methods Based on Operator Type 

For the jump models considered in Section 2.3 the semidiscrete matrix A can be 
written in the form 

A = D + J, (23) 

where D and J correspond to the differential operator and integral operator, respec¬ 
tively. The matrix D is sparse while in general J is a dense matrix or has dense 
blocks. In view of the different nature of these two matrices it can be preferable to 
employ an operator splitting method based on them. 

In [3], Andersen and Andreasen describe a generalized 0-method 
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(I - dDAt\y - J) t/„ = (I + (1 - 0d) 4fD + (1 - dj) Ati) Un- 1 (24) 

assuming here G = 0. The standard choice Qo = I and = 0 corresponds to the 
IMEX Euler method: it treats the stiff differential part implicitly, using the back¬ 
ward Euler method, and the nonstiff integral part explicitly, using the forward Euler 
method. This choice yields first-order consistency. The benefit is that it is not neces¬ 
sary to solve dense linear systems involving the matrix J. Instead, in each time step 
only one multiplication with J is required. This approach has been considered and 
analysed in [17]. 

In [26] an extrapolation approach is advocated based on the IMEX Euler method. 
Here approximations at a given fixed time are computed for a decreasing sequence 
of step sizes and then linearly combined so as to achieve a high order of accuracy. 

In [3] second-order consistency is obtained through an alternating treatment of 
the D and J parts. They propose to take a At jl substep with = 1 and 6j = 0 
followed by a At/2 substep with 9d = 0 and 6j = 1. Here linear systems involving 
the dense matrix J need to be solved, for which the authors employ EET 

In [22] the original 0-method is analyzed, where the linear system in each time 
step is solved by applying a fixed-point iteration on the jump part following an idea 
in [77]. 

The following, second-order/MEX midpoint scheme has been considered in e.g. 
[26, 57, 58, 72], 

{l-AtD)Un = (I-f 4fD)t4_2 + 24fJt/„_i-f24fG„_i. (25) 

The scheme (25) can be viewed as obtained from the semidiscrete system (17) at 
f„_i by the approximations Dt/„_i Ri jiy{Un -b t/„_ 2 ) and t/„_i ^ Un-i). 

Two subsequent second-order IMEX methods are the IMEX-CNAB scheme 

(l-fD)t/„ = (l-ffD)t/„_i + fj(3t/„_i-t/„_2)+4fG„_i/2 (26) 

and the IMEX-BDF2 scheme 

(|I - AtD) Un = 2Un-X - \Un-2 + Ati (2t/„_i - Un-l) + AtCn- (27) 

These schemes have recently been applied for option pricing in [73] and can be 
regarded as obtained by approximating the semidiscrete system (17) at f„_i /2 = 
5 (f„ + fn- 1 ) and at f„, respectively. 

The IMEX schemes (25), (26), and (27) were studied in a general framework, 
without application to option valuation, in [28]. Here it was noted that such schemes 
can be considered as starting with an implicit method and then replacing the nonstiff 
part of the implicit term by an explicit formula using extrapolation based on previous 
time steps. An overview of IMEX methods is given in [48]. 

In general, IMEX methods are only conditionally stable, that is, they are stable 
for a sufficiently small time step At. Eor example, the IMEX midpoint scheme (25) 
and the IMEX-CNAB scheme (26) are stable whenever XAt < 1 and the Xu term 
in (10) is included in D; see [73]. Recall that X denotes the jump activity. 
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The schemes discussed in this section are of the linear multistep type. For IMEX 
schemes of Runge-Kutta type applied to jump models we mention [10]. 


5.4 Operator Splitting Method for Linear Complementarity 
Problems 

The fully discrete LCPs obtained by spatial and temporal discretization of (14) for 
American-style options are more difficult to solve than the corresponding systems 
of linear equations for the European-style counterparts. It is desirable to split these 
LCPs into simpler subproblems. Here we describe the operator splitting method 
considered in [49, 53] which was motivated by splitting methods for incompressible 
flows [13, 31]. To this purpose, we reformulate LCPs with Lagrange multipliers. 

The 0-method discretization (18) naturally gives rise to the following, fully dis¬ 
crete LCP 


Bt/„-Ct/„_i-AfG„_i+0>O, 

Un > t/o, (Bt/„ - CUn-l - AtGn-x+ef {U„ - Uo) = 0, 


(28) 


where B = I — 9 At A, C = I -F (1 — 0 )AtA, and A is assumed to be constant in time. 
By introducing a Lagrange multiplier vector X,,, the LCP (28) takes the equivalent 
form 


Bt/„ - Ct/„_i - AfG„_i+e = AtX„ > 0, 
Un>Uo, {Xnf{Un-Uo)=0. 


(29) 


The basic idea of the operator splitting method proposed in [49] is to decouple in 
(29) the first line from the second line. This is accomplished by approximating the 
Lagrange multiplier A„ in the first line by the previous Lagrange multiplier A„_ i. 
This leads to the system of linear equations 


Bt7„ = CUn-i -f-AfG„_i+0 +AtX„-i. 


(30) 


After solving this system, the intermediate solution vector Un and the Lagrange 
multiplier A„ are updated to satisfy the (spatially decoupled) equation and comple¬ 
mentarity conditions 


( Un-Un = At[Xn-X„^x), 

{ T (31) 

[A„>0, G„>Go, (A„)^ (G„ - Go) = 0. 

Thus, this operator splitting method for American options leads to the solution of 
linear systems (30), which are essentially the same as for European options, and a 
simple update step (31). This update can be performed very fast, at each spatial grid 
point independently, with the formula 
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{ i^nj 1,1 ; Oj , if Ufj^i AtXf^—i i > Uqj , 

; / ~ 

{Uo,i , Xn-ij + ^ [Uo^i - u„jj j , otherwise. 

The above operator splitting approach has been studied for more advanced time 
discretization schemes of both linear multistep and Runge-Kutta type in [49, 53]. 
Moreover, it has recently been effectively combined with IMEX schemes in [72] for 
the case of jump models and with ADI schemes in [37] for the case of the Heston 
model. For instance, the pertinent adaptations of the IMEX-CNAB scheme and the 
MCS scheme are 

(l— 4^D) Un = (l + 4^D) Un-\ + — Un-l) + AtG^_ij2 + At X,i-\, 

and 

' To = t/„_i + + Af 

Yj = + dAt {Yj{tn,Yj) - (7=1,2,...,^), 

Yo = Yo + 9 At (Fo(f„,Fj.) — Fo(f„_i,t/„_i)), 

< 

fb = To + (i - e)At iF{tn,Yk)-F{t„-i,U„-i)), 

Yj = F,_i + 9 At (Fjitnjj) - Fj{t„-uU„-i)) (7=1,2,...,^), 

, U„ = 

respectively, followed by the update (32). The other three ADI schemes from Sec¬ 
tion 5.2 are adapted analogously. Note that only a Af i term has been added to the 
first line of the MCS scheme (21). Accordingly, like for the 0-method, the amount 
of computational work per time step is essentially the same as for the corresponding 
European-style option. 


6 Solvers for Algebraic Systems 

The implicit time discretizations described in Section 5 lead, in each time step, to 
systems of linear equations of the form 




(33) 


or LCPs of the form 


Bt/ > •?, U><t>, 

{m-'Ff {u-<p) = o 


( 34 ) 
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with given matrix B and given vectors <P, 'P. For models without jumps, semidis¬ 
cretization by finite difference, finite volume, and finite element methods yields 
sparse matrices B. For one-dimensional models, the central FDs (15) and (16) lead 
to tridiagonal B. For higher dimensional models they give rise to matrices B with a 
large bandwidth whenever classical (non-splitted) time stepping schemes are ap¬ 
plied. On the other hand, for the operator splitting methods based on direction 
(cf. Section 5.2) one also acquires tridiagonal matrices (possibly after renumber¬ 
ing the unknowns). Wider FD stencils lead to additional nonzero diagonals. Time 
discretization of jump models with an implicit treatment of jumps makes B dense. 


6.1 Direct Methods 

The system of linear equations (33) can be solved by a direct method using LU 
decomposition. This method first forms a lower triangular matrix L and an upper 
triangular matrix U such that B = LU. After this the solution vector U is obtained 
by solving first L,V = P and then Ut/ = V. 

Let m denote the dimension of the matrix B. For tridiagonal B, or more generally 
matrices with a fixed small bandwidth, the LU decomposition yields optimal com¬ 
putational cost in the sense that the number of floating point operations is of order m. 
Flence, it is very efficient for one-dimensional models and for higher-dimensional 
models when operator splitting schemes based on direction are applied. 

For two-dimensional models with classical time stepping schemes, a LU decom¬ 
position can be formed by order floating point operations if a nested dissection 
method can be used and then the computational cost of the solution is of order 
mlogm, see [21, 29]. For higher-dimensional models with classical time stepping 
schemes, the computational cost is less favorable. 

For a general matrix B, solving the LCP (34) requires iterative methods. How¬ 
ever, in the special case that B is tridiagonal, the solution vector satisfies t/, = <Pi 
(1 < / < ;o), Ui > <l>i (i'o < i < m) for certain /q and some additional assumptions 
hold, the Brennan-Schwartz algorithm [9] gives a direct method to solve the LCP; 
see also [1,51, 55]. After inverting the numbering of the unknowns to be from right 
to left, represented by a permutation matrix P, this algorithm is equivalent to apply¬ 
ing the LU decomposition method to the corresponding linear system with matrix 
PBP where the projection step is carried out directly after computing each compo¬ 
nent in the back substitution step with U. More precisely the back substitution step 
reads after the renumbering of unknowns: 

f tJm — IttaxIVm/Um,/?! 7 7 

< (35) 

[ Ui = max{(y, - U,-,+it/ 7 +i)/U,v ,0,} (/ = m - 1 ,m - 2 ,..., 1). 

The Brennan-Schwartz algorithm is essentially as fast as the LU decomposition 
method for linear systems and, thus, it has optimal computational cost. 
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6.2 Iterative Methods 

There are many iterative methods for solving systems of linear equations. The 
two most important method categories are the stationary iterative methods and the 
Krylov subspace methods. Well-known Krylov subspace methods for the, typically 
unsymmetric, system (33) are the generalized minimal residual (GMRES) method 
[70] and the BiCGSTAB method [84]. In the following we discuss a stationary it¬ 
erative method in some more detail which is familiar in finance applications. The 
successive over-relaxation (SOR) method reads 



(36) 


for i = 1,2,... ,m, k = 0,1,2,..., where o is a relaxation parameter. This method 
reduces to the Gauss-Seidel method in the case (0=1. The convergence rate of 
the iteration (36) can be improved significantly by an optimal choice of (O. Still the 
number of iterations to reach a given accuracy typically grows with m, that is, when 
the spatial grid is refined the convergence slows down. 

The SOR iteration can be generalized for LCPs by performing a projection after 
each update [19]; see also [30]. This method is called the projected SOR (PSOR) 
method and it reads 



(37) 


(/ = 1,2,..., OT, k = 0,1,2,...). As can be expected, the PSOR method suffers from 
the same drawback as the SOR method mentioned above. 


6.3 Multigrid Methods 

The aim of multigrid methods for solving linear systems (33) is to render the number 
of iterations essentially independent of the problem size m. The stationary iterative 
methods typically reduce high frequency errors quickly, while low frequency errors 
are reduced much more slowly. The idea of multigrid methods is to compute effi¬ 
ciently corrections to these slowly varying errors on coarser spatial grids. The multi¬ 
grid methods can be divided into geometrical and algebraic methods. With the ge¬ 
ometrical methods discretizations are explicitly constructed on a sequence of grids 
and transfer operators between these grids are explicitly defined. Algebraic multi¬ 
grid (AMG) methods [69, 76] build the coarse problems and the transfer operators 
automatically using the properties of the matrix B. The details of these methods are 
beyond the scope of this chapter and we refer to e.g. [82] for details and extensive 
literature on this. 
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Several versions of multigrid methods also exist for LCPs. Brandt and Cryer in¬ 
troduced in [8] a projected full approximation scheme (PFAS) multigrid method 
for LCPs. American options under stochastic volatility were priced using the PFAS 
method in [15, 65]. A projected multigrid (PMG) method for LCPs introduced in 
[68] resembles more closely a classical multigrid method for linear problems. This 
method has been used to price American options in [52, 68]. Recently, an AMG 
method was generalized for LCPs in [81]. The resulting method is called the pro¬ 
jected algebraic multigrid (PAMG) method and resembles the PMG method in the 
treatment of the complementarity conditions. 


7 Numerical Illustrations 

In the following we price European and American put options under a hierarchy of 
models: Black-Scholes, Merton, Heston, and Bates. The interest rate, the maturity 
time, and the strike price are always taken as 

r = 0.03, T = Q.5, and /ST =100. 

For the purpose of illustration. Fig. 1 and Fig. 2 show fair values of European and 
American options, respectively, under the four considered models with the model 
parameters described in the following sections. 



Fig. 1 The fair values of European put options for the asset prices 75 < i < 125 and the volatility 
O' = 0.2 (the variance v = 0.04) under the four considered models. 
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Fig. 2 The fair values of American put options for the asset prices T5 < s < 125 and the volatility 
a = 0.2 (the variance v = 0.04) under the four considered models. 


7.1 Black-Scholes model 

In the case of the Black-Scholes model, we price American put options. The volatil¬ 
ity in the model (1) is taken as 

(7 = 0.2 

and the following boundary conditions are employed; 

u{0,t) = K for 0<t<T, (38) 

Us{Smax,t)=0 for Q<t<T. (39) 

The Neumann boundary condition (39) introduces a modeling error as it is not ex¬ 
actly fulfilled by the actual option price function. If S'max is taken sufficiently large, 
however, this error will be small in the region of interest. 

For the spatial discretization of the Black-Scholes PDF (2), we apply FD formu¬ 
las on nonuniform grids such that a large fraction of the grid points lie in the region 
of interest, that is, in the neighborhood of 7 = /T. 

For the construction of the spatial grid we adopt [36]. Let integer m i > 1, constant 
c> 0, and 0 < 5ieft <K < ^ught < be given. Let equidistant points < 

< ... < = ^max be given with distance and 
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^min = sinh ^ : 

^ _ ‘bright ‘^left 

Sint ~ 1 

C 

e e , • I - 1 / ‘^max ~ ‘bright X 

Smax = Sint + Smh ( ---^ 1 . 

Then we define a nonuniform grid 0 = sq < ^i < ■■ ■ < Sm^ = 5max by the transfor¬ 
mation 

Si = (p{^i) (0</<mi), (40) 

where 

pieft +c-sinh(^) min — ^ 0 ), 

<P(^)= < ^left +C-^ (0<<^ <<^int), 

‘bright “I" ^ ‘ '^int) (‘Dint — ‘Dmax)- 

The grid (40) is uniform inside [5ieft,5right] and nonuniform outside. The parameter 
c controls the fraction of grid points i, that lie inside [5ieft,5right]. The grid is smooth 
in the sense that there exist real constants Co,Ci,C 2 > 0 such that the grid sizes 
Asi = Si — s,_i satisfy 

CoA^ <Asi<CiA^ and |Ai,+i - Z\i,| < C 2 (4^)^ (41) 

uniformly in i and nii. For the parameters in the grid we make the (heuristic) choice 

5max = 8 X, c = ^, 5ieft = , bright = min X. 

The semidiscretization of the initial-boundary value problem for the Black- 
Scholes PDF is then performed as follows. At the interior grid points each spatial 
derivative appearing in (2) is replaced by its corresponding second-order central FD 
formula described in Section 4. At the boundary s = S'max the Neumann condition 
(39) gives dufds. Next, d^ujds^ is approximated by the central formula with the 
value at the virtual point 5max + Asmi defined by linear extrapolation using (39). 

Concerning the initial condition, we always replace the value of the payoff func¬ 
tion 0 at the grid point s, nearest to the strike K by its cell average, 

1 /’■'■/+ 1/2 

- / max(X — 
h 

where 

■Si- 1/2 = +■5/); ■S(+l /2 = 5 ('S!+'Si+l)i h = ^i+ 1/2 —^i- 1 / 2 - 

This reduces the dependency of the discretization error on the location of the strike 
relative to the s-grid, see e.g. [77]. 





22 


Karel in’t Hout and Jari Toivanen 


The time discretization is performed by the Crank-Nicolson method with Ran- 
nacher smoothing. The time stepping is started by taking two backward Euler steps 
using the time step \At. With this choice all time steps are performed with the 
same coefficient matrix I — jAtA. Furthermore, halving the time step with the Eu¬ 
ler method helps to reduce the additional error caused by this method. Note that we 
count these two Euler steps as one time step in order to keep the notations conve¬ 
nient. 

We define the temporal discretization error to be 

E(mi,iV) = max { \UN,i - Ui{T)\ : {K < st < ^K} , (42) 

where Unj denotes the component of the vector Un associated to the grid point s,. 
We study the temporal discretization errors on the grids {m\,N) = (160,2^) for 
k = 0,1,..., 10. The reference price vector U{T) is computed using the space- 
time grid (160,5000). Fig. 3 compares the temporal errors of the smoothed Crank- 
Nicolson method with and without the operator splitting method for LCPs described 
in Section 5.4. For larger time steps the Crank-Nicolson method without splitting is 
more accurate. In this example the convergence rate of the splitted method is slightly 
less than second-order and a bit higher than the convergence rate of the unsplitted 
method. Thus, for smaller time steps the operator splitting method is slightly more 
accurate. 



Fig. 3 The temporal discretization errors for the American option under the Black-Scholes model 
for the smoothed Crank-Nicolson method with and without the operator splitting method for LCPs. 
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7,2 Merton model 

Under the Merton jump diffusion model, we price European and American put op¬ 
tions. For the jump part of the model, the jump activity, the mean of the normal 
distribution, and its standard deviation are taken as 


1=0.2, 5 = 0.4, and 7 =-0.5 


(43) 


respectively; see ( 8 ). The boundary condition at i = 0 is given by (5) for the Euro¬ 
pean put option and by (38) for the American put option. At the truncation boundary 
s = 5max, we use the Neumann boundary condition (39). 

The same space-time grids are considered as with the Black-Scholes model in 
Section 7.1 and also the spatial derivatives are discretized in the same way. For the 
integral term, we use a linear interpolation for u between grid points and take u to 
be zero for s > Smax- The formulas for the resulting matrix J are given in [71], for 
example. 

For the time discretization, we apply the IMEX-CNAB scheme, which is always 
smoothed by two Euler steps with the time step jAt. In these first steps the backward 
Euler method is used for the discretized differential part D and the forward Euler 
method is used for the discretized integral part J. For European options, these steps 
are given by 


(I - f D) t/ 1/2 = t/o + f Jt/o + f Gi/ 2 , 

(I - f D) Ui = t/i /2 + /2 + f Gi. 


In the absence of jumps, these steps reduce to the same Rannacher smoothing used 
with the Black-Scholes model. After these two steps the IMEX-CNAB scheme 
defined by (26) is employed. 

We study the temporal discretization errors for European and American options 
on the same grids (mi,A) = (160,2^), k = 0 , 1 ,..., 10 , and using the same error 
measure (42) as before. Fig. 4 shows the temporal errors for the European option us¬ 
ing the IMEX-CNAB scheme and the Crank-Nicolson method with classical Ran¬ 
nacher smoothing. We observe that the temporal errors for the two methods are 
essentially the same and they exhibit second-order convergence. 

Fig. 6 shows the same temporal errors for American options using the IMEX- 
CNAB scheme with operator splitting for LCPs and the Crank-Nicolson method 
without splitting. The convergence result for the two methods is very similar to the 
case of the Black-Scholes model in Section 7.1. Thus, for larger time steps the 
Crank-Nicolson method is more accurate while for smaller time steps the IMEX- 
CNAB scheme with splitting is more accurate. 

In order to gauge the effectiveness of the proposed discretizations, we report the 
total discretization errors for the European option on the space-time refining grids 
(mi, A) = 2^(10,2), k = 0,1,..., 6 . The total discretization error is defined by 


e(mi,A) = max{ \UN,i — u{si,T) \ : \K < Si < . 


(44) 
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The reference price function ii is computed on the space-time grid (10240,2048). 
Fig. 5 shows the total error for the European option using the IMEX-CNAB scheme 
and the Crank-Nicolson method. As with the temporal errors the total errors for 
both methods are essentially the same and both show a second-order convergence 
behavior. 



Fig. 4 The temporal discretization errors for the European option under the Merton model with 
the IMEX-CNAB scheme and the Crank-Nicolson method, both with smoothing. 


7.3 Heston model 

Under the Heston stochastic volatility model we consider European and American 
put options as well. Eor the mean-reversion rate, the long-term mean, the volatility- 


of-variance and the correlation the following values are taken: 

k = 2, 77 =0.04, (7 = 0.25, and p = —0.5. (45) 

The spatial domain is truncated to [0,5max] x [0, Vmax] with 5max = and Vmax = 5. 
The following boundary conditions are imposed: 

u{0,v,t) = df ■ K for 0 < V < Umax, 0 < f < r, (46) 

Mi(5max, V, f) = 0 for 0 < V < Umax, 0 <. t ^ T, (47) 

Mv(7,Vmax,f) = 0 for 0 < S < 5niax, 0 < f < T, (48) 
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Fig. 5 The total discretization errors for the European option under the Merton model with the 
IMEX-CNAB scheme and the Crank-Nicolson method, both with smoothing. 



Fig. 6 The temporal discretization errors for the American option under the Merton model with 
the IMEX-CNAB scheme together with the operator splitting method for LCPs, and the Crank- 
Nicolson method, both with smoothing. 
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where df = in the European case and df = 1 in the American case. At the 
degenerate boundary v = 0 the Heston PDE holds in the European case and it is 
assumed that the Heston LCP holds in the American case. The two conditions at 
■s = 5'max and v = introduce a modeling error, as they are not exactly fulfilled 
by the actual option price function, but in our experiments this error is small on the 
region of interest in the (s, v)-domain. 

Eor the spatial discretization of the Heston PDE and Heston LCP we apply ED 
formulas on Cartesian grids. Here nonuniform grids are used in both the s- and v- 
directions such that a large fraction of the grid points lie in the neighborhoods of 
s = K and v' = 0, respectively. This is the region in the (i,v')-domain where one 
wishes to obtain option prices. Next, the application of such nonuniform grids can 
greatly improve the accuracy of the ED discretization as compared to using uniform 
grids. This is related to the facts that the initial function (4) possesses a discontinuity 
in its first derivative at s = /T and that for v « 0 the Heston PDE is convection- 
dominated. The grid in the s-direction is taken identical to that in Section 7.1. 

To construct the grid in the v-direction, let integer m2 > 1 and constant d >0 and 
let equidistant points be given by '^fj = j ■ for 7 = 0,1 ,..., m2 with 



Then a smooth, nonuniform grid 0 = v'o < vq < ... < v„i^ = Vmax is defined by 


(f-sinh(v/7) {0<j<m2). 


(49) 


The parameter d controls the fraction of grid points Vj that lie near v = 0. We heuris- 
tically choose 



The semidiscretization of the initial-boundary value problem for the Heston PDE 
and Heston LCP is performed as follows. In view of the Dirichlet condition (46), the 
grid in [0,5max] x [0,Vmax] is given by {{si,Vj) : 1 < / < mi, 0 < 7 < m2}. At this 
grid, each spatial derivative is replaced by its corresponding second-order central 
ED formula described in Section 4 with a modification for the boundaries v = 0, 
^ — ‘^max5 and V = Eniax- 

At the boundary v = 0 the derivative du/dv is approximated using a second-order 
forward formula. All other derivative terms in the v-direction vanish at v = 0, and 
therefore do not require further treatment. 

At the boundary s = 5max the spatial derivatives in the i-direction are dealt with 
as in Section 7.1. Note that the Neumann condition (47) at i = 5max implies that the 
mixed derivative d^ujdsdv vanishes there. 

At the boundary v = Vmax the spatial derivatives in the v-direction need to be 
considered. This is done fully analogously to those in the s-direction at s = Smax 
using now the Neumann condition (48). 

Define the temporal discretization error by 
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= max{ — t//(7’)| : jK < Sj < jK, 0 < vj < 1} , (50) 

where the index I corresponds to the grid point The reference vector 

U{T) is computed using {mi,m 2 ,N) = (160,80,5000). We study these errors for 
{m\,m2,N) = (160,80,2^) with A: = 0,1,..., 10 and three methods: the Do scheme 
with 0 = 5 and smoothing, the MCS scheme with 0 = 5 without smoothing, and 
the Crank-Nicolson scheme with smoothing. 

Fig. 7 displays the obtained results for the European put option. As a first ob¬ 
servation, for all three methods the temporal errors are bounded from above by a 
moderate value and decrease monotonically as N increases. The error graphs for the 
MCS and Crank-Nicolson schemes are almost identical and reveal a second-order 
convergence behavior. The Do scheme only shows first-order convergence. Clearly, 
the convergence orders observed for the three methods agree with their respective 
classical orders of consistency. Additional experiments by substantially changing 
(mi,m 2 ) indicate that for all three methods the temporal errors are almost unaf¬ 
fected, which is a desirable property and suggests convergence in the so-called stijf 
sense. Whereas their results are not displayed, we mention that the CS scheme with 
0 = 5 and smoothing and the HV scheme with 0 = 5 - 1 - gv/S without smoothing 
behave similarly to the MCS scheme in this experiment, with slightly larger errors. 

Fig. 9 displays the obtained results for the American put option. Our observa¬ 
tions are analogous to those made above in the case of the European option. It is 
interesting to note, however, that the Do scheme often has temporal errors that are 
almost the same as for the MCS and Crank-Nicolson schemes. But if N gets suffi¬ 
ciently large, then a first-order convergence behavior for this method indeed sets in. 
For the Crank-Nicolson scheme a small deviation from second-order is seen when 
N is large. This disappears however when other values (mi, m 2 ) are considered. Ad¬ 
ditional experiments by substantially changing (mi,m 2 ) indicate that for all three 
methods the temporal errors are at most mildly affected. 

We next consider, in the European put option case, the total discretization error 
defined by 

e{m\,m2,N) =max{\UNj — u{si,Vj,T)\ : < Si < jK, 0 < vj < l} , (51) 

with index I corresponding to the grid point {si,Vj). Here exact solution values u 
are computed by a suitable implementation of Heston’s semi-closed form analytical 
formula [38]. Note that the modeling error, which is due to the truncation of the 
domain of the Heston PDF to a bounded set, is also contained in e{mi,m2,N). In 
our experiment, this contribution is negligible. 

Fig. 8 displays the total discretization errors for (mi,m 2 ,N) = 2^(10,5,2) with 
k = 0,1,... ,6 and the three schemes under consideration in this section. With the 
MCS and Crank-Nicolson schemes the total errors are essentially the same and a 
second-order convergence behavior is observed. With the Do scheme, the total errors 
are almost same as these two schemes up to k = 4, but then the convergence drops 
to the expected first-order. 
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For a more extensive numerical study of ADI schemes in the (two-dimensional) 
Heston model we refer to [39] for European-style options and to [37] for American- 
style options. For three-dimensional PDEs in finance, such as the HHW PDF, the 
numerical convergence of ADI schemes has been investigated in [35, 36] and for 
a four-dimensional PDE in [34]. In these references a variety of parameter sets has 
been considered, including long maturity times and cases where the Feller condition 
is strongly violated, together with various barrier options and the approximation of 
hedging quantities. 



Fig. 7 Temporal discretization errors in the case of the European put option under the Heston 
model. The time discretization methods are: the Do scheme with 9 = 5 and smoothing, the MCS 
scheme with 9 = 5 without smoothing, and the Crank-Nicolson scheme with smoothing. 


7.4 Bates model 

We price European and American put options under the Bates model. The boundary 
conditions are given by (46)-(48). For the stochastic volatility part of the model the 
parameters are taken the same as for the Heston model and they are given by (45). 
For the jump part, the parameters are the same as for the Merton model and they 
are given by (43). The discretizations are based on the same grids and the spatial 
derivatives are discretized in the same way as with the Heston model in Section 
7.3. For the jump integral, the same discretization is used as with the Merton model 
in Section 7.2. We consider here the IMEX-CNAB scheme and Crank-Nicolson 
method both applied with smoothing as for the Merton model. 
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Fig. 8 Total discretization errors in the case of the European put option under the Heston model. 
The time discretization methods are: the Do scheme with 6 = 5 and smoothing, the MCS scheme 
with 6 = 5 without smoothing, and the Crank-Nicolson scheme with smoothing. 



Fig. 9 Temporal discretization errors in the case of the American put option under the Heston 
model. The time discretization methods are: the Do scheme with 6 = 5 and smoothing, the MCS 
scheme with 9 = 5 without smoothing, and the Crank-Nicolson scheme with smoothing. 
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As with the Heston model, we consider the temporal discretization errors on the 
grids (nil,m 2 ,N) = (160,80,2^), A: = 0,1,..., 10. The reference price vector U(T) 
is computed using the space-time grid (160,80,5000). The temporal discretization 
errors e(mi,m 2 ,N) are shown for the European option in Fig. 10 and for the Amer¬ 
ican option in Fig. 12. The plots show the errors for the IMEX-CNAB scheme and 
the Crank-Nicolson method. For the American option the operator splitting method 
for FCPs is used with the IMEX-CNAB scheme. As with other models, the tempo¬ 
ral errors for the European option are very similar for both methods and they both 
exhibit second-order convergence. For the American option, the difference between 
the methods is less pronounced than with the Black-Scholes and Merton models. 
Still the Crank-Nicolson method is slightly more accurate than the operator split¬ 
ting method for large time steps and the reverse is true for small time steps. In this 
example the convergence rates seem to be between 1.5 and 2.0. 

We computed the total discretization errors e(mi,m 2 ,N) for the European option 
on the grids {mi,m 2 ,N) = 2^(10,5,2), k = 0,1,... ,6. The reference prices are com¬ 
puted on the space-time grid (2560,1280,512). Fig. 11 shows the total errors for the 
IMEX-CNAB scheme and the Crank-Nicolson method. As with the other models, 
the total errors for both methods are virtually the same and both have second-order 
convergence of the total error. 



Fig. 10 The temporal discretization errors for the European option under the Bates model with the 
IMEX-CNAB scheme and the Crank-Nicolson method, both with smoothing. 
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Fig. 11 The total discretization errors for the European option under the Bates model with the 
IMEX-CNAB scheme and the Crank-Nicolson method, both with smoothing. 



Fig. 12 The temporal discretization errors for the American option under the Bates model with 
the IMEX-CNAB scheme together with the operator splitting method for LCPs, and the Crank- 
Nicolson method, both with smoothing. 
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8 Conclusions 

We have discussed numerical solution methods for financial option valuation prob¬ 
lems in the contemporary partial(-integro) differential equations framework. These 
problems are often multidimensional and can involve nonlocal integral operators 
due to jumps incorporated in the underlying asset price models. The early exercise 
feature of American-style options gives rise to linear complementarity problems, 
which are nonlinear. All these properties add complexity to the discrete problems 
obtained by classical implicit numerical methods and renders their efficient solution 
a challenging task. The efficient computation of option values is, however, crucial 
in many applications. In this chapter an overview has been given of various types 
of operator splitting methods for the discretization in time, which yield in each time 
step a sequence of discrete subproblems that can be handled much more easily and 
efficiently without essentially influencing the accuracy of the underlying discretiza¬ 
tion. The following highlights the different operator splitting methods presented in 
this chapter. 

For multidimensional models the directional splitting methods of the ADI type 
offer a fast, accurate, and easy-to-implement way for the numerical time stepping. 
They are adapted to effectively deal with mixed spatial derivative terms, which are 
ubiquitous in finance. ADI schemes lead to a sequence of sparse linear subproblems 
that can be solved by LU decomposition at optimal computational cost, that is, the 
number of required operations is directly proportional to the number of unknowns. 
The MCS and HV schemes, with a proper choice of their parameter 0, are recom¬ 
mended as these show stability and second-order convergence and reveal a better 
inherent smoothing than second-order CS. 

The spatial discretization of jumps models for the underlying asset price yields 
dense matrices. All classical implicit time discretization schemes require solving 
systems with these dense matrices. By employing an IMEX method like the IMEX- 
CNAB scheme advocate here, with an explicit treatment of (finite activity) jumps 
and an implicit treatment of the remainder of the operator, each time step involves 
only multiplications with these dense matrices. This is computationally a much eas¬ 
ier task and can be often performed very fast using EET The accuracy and stability 
of the IMEX-CNAB scheme are good when the jump activity is not very high, e.g. 
less than several jumps per year. 

Iterative methods like the PSOR method for solving LCPs resulting from the 
pricing of American-style options often converge slowly. We discussed an operator 
splitting method based on a Lagrange multiplier formulation, treating in each time 
step the early exercise constraint and complementarity condition in separate sub¬ 
problems, where the main subproblem is essentially the same as for the European- 
style counterpart. With this approach it is easy to adapt a European option pricer to 
American options. We presented such an adaptation for ADI and IMEX methods. 
Also, it is applicable for most models of underlying asset prices. Numerical experi¬ 
ence with this operator splitting method indicates that the accuracy stays essentially 
the same as in the case of the original LCP, but there can be a major reduction in 
computational time. 
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